function Z = geval(dof_map,V,T,d,c,X,Y)
% This function evaluates the piecewise poly over the triangulation [V,T] defined 
% by the vector c of Bnet coeffs at the points with coordinates X and Y,
% we'd better make it a quick version or translate it into c or fortran.
Z = nan*ones(size(X));
for tri = 1:size(T,1)
    [b1,b2,b3] = bary(V(T(tri,1),:),V(T(tri,2),:),V(T(tri,3),:),X(:),Y(:));
    I = find(b1>=0 & b2>=0 & b3>=0);
    if ~isempty(I)
        Z(I) = vdm23(d,b1(I),b2(I),b3(I))*c(dof_map(:,tri));
    end;
end
end

function [lam1,lam2,lam3] = bary(V1,V2,V3,X,Y)
One = ones(size(X(:)'));
A = [1 1 1 ;V1(1),V2(1),V3(1);V1(2),V2(2),V3(2)];
lam = A\[One;X(:)';Y(:)'];
lam1 = reshape(lam(1,:),size(X));
lam2 = reshape(lam(2,:),size(X));
lam3 = reshape(lam(3,:),size(X));
end


function Mat = vdm23(d,b1,b2,b3)
% comput the Bform of degreed d at points (b1,b2,b3);
m = (d+1)*(d+2)/2;
plot_m = length(b1);
[I,J,K] = indices(d);
IM = diag(I)*ones(m,plot_m);
JM = diag(J)*ones(m,plot_m);
KM = diag(K)*ones(m,plot_m);

plot_IM = diag(b1)*ones(plot_m,m);
plot_JM = diag(b2)*ones(plot_m,m);
plot_KM = diag(b3)*ones(plot_m,m);
Mat = (plot_IM).^(IM').*(plot_JM).^(JM').*(plot_KM).^(KM');
IF = gamma(I+1);
JF = gamma(J+1);
KF = gamma(K+1);
A = factorial(d)*ones(plot_m,m)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;
end